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An explicit optimal linear spatial predictor is derived. The spatial correlations are imposed by 
means of Gibbs energy functionals with explicit coupling coefficients instead of covariance matrices. 
The model inference process is based on physically identifiable constraints corresponding to distinct 
terms of the energy functional. The proposed predictor is compared with the geostatistical linear 
optimal filter (kriging) using simulated data. The agreement between the two methods is excellent. 
The proposed framework allows a unified approach to the problems of parameter inference, spatial 
prediction and simulation of spatial random fields. 
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INTRODUCTION 

Spatial prediction of physical variables from samples 
that are irregularly distributed in space is a task with 
applications in many fields of science and engineering, 
including subsurface hydrology 0, 01 j oil reservoir engi- 
neering |3, yh environmental pollutant mapping and risk 
assessment 5], mining exploration and reserves estima- 
tion , environmental health studies Q , image analysis 
and neuroscience Physical quantities of econom- 
ical and environmental interest include mineral grades, 
concentrations of environmental pollutants, soil and rock 
permeability and flow fields in oil reservoirs. Modeling 
the variability of such processes is based on the theory 
of spatial random fields (SRFs) |lf|. Knowledge of spa- 
tial correlations in SRFs enables (i) generating predic- 
tive isolevel maps (ii) estimating prediction uncertainty 
and (iii) developing simulations that reconstruct proba- 
ble scenarios conditioned on the data. The classical ap- 
proach is based on Gaussian SRF's (GSRF's) and various 
generalizations for non-Gaussian distributions [ill fl2| . 
For GSRF's the spatial structure is determined from the 
covariance matrix, which is estimated from the available 
sample (s). 

Let fi £ M. d denote the area of interest, and |0| its vol- 
ume. An SRF state (realization) in Q can be decomposed 
into a deterministic trend m x (s), a correlated fluctua- 
tion X\(s), and an independent random noise term, e(s), 
i.e., X(s) — m x (s) + X\(s) + e(s). The trend represents 
large-scale variations obtained in principle by ensemble 
averaging, i.e. m x (s) = E[X(s)]. In practice, the trend 
is often determined from a single available realization. 
The fluctuation represents 'fast variations' describing fine 
structure above the resolution limit A. The random noise 
represents non-resolved inherent variability, purely ran- 
dom additive noise, or non-systematic measurement er- 
rors. The fluctuation typically is a second-order station- 
ary SRF, or an intrinsic SRF with second-order station- 



ary increments [loj . The residual SRF after trend re- 
moval is a zero-mean fluctuation: X*(s) = A^a(s) + e(s). 
The inference process focuses on determining a model for 
ATa(s) from a set of possibly noisy observations X*(s). 

In statistical physics the probability density function 
(pdf) of any fluctuation field X(s) governed by an en- 
ergy functional H[X(s)] is expressed as f x [X(s)] = 
Z~ x exp {— if[X(s)]} , where Z is the partition function. 
In classical geostatistics, the Gaussian joint pdf for a set 
of fluctuations X = {X(si),i — 1 . . . , N} is expressed in 
terms of H[X] = ±X( Si ) [C^ 1 X( Sj ), where [CJr 1 is 
the inverse covariance matrix, and summation is implied 
over repeated indices. Instead, Spartan Spatial Random 
Fields (SSRF's) model spatial correlations in terms 
of 'interactions'. This change in viewpoint has important 
consequences for model inference and spatial prediction. 



THE FGC MODEL 

In the fluctuation-gradient-curvature (FGC) SSRF 
model is defined and its properties investigated. The 
continuum FGC model involves the following if[X]: 



Hf gc [X\ 



2r h 



^djds [S (b) + vi e Si(s) + e S 2 (s)] , 

(1) 

where S (s) = [X x (s)} 2 , 5i(s) = [VX A (s)] 2 , and S 2 (s) = 
[V 2 X\(s)Y . The model is characterized by four param- 
eters: the scale coefficient rjo, the covariance-shape co- 
efficient the characteristic length £, and the cutoff 
wavevector fc max . Bochner's permissibility theorem || 
for the positive definiteness of the covariance function 
requires r\\ > —2 if fc max — > oo. In statistical physics 
terminology, l x = Vl ^ d /^Vo) and £ 2 = ^ d /^Vo) 
represent the coupling strengths of the gradient and cur- 
vature terms. A coarse-graining kernel is used to cut 
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off the fluctuations at /c max 0, 0], leading to band- 
limited covariance spectral density. If k max is finite, the 
field's configurations are almost surely differentiable [l4| . 
If fc max is infinite, generalized gradient and constraints 
should be used. The coarse-graining kernel implies that 
the SRF X x is a generalized SRF fig . 

A moment-based method for parameter estimation was 
proposed and validated with simulated data The 
inference process is based on matching ensemble con- 
straints E[Sj(s)] with their sample counterparts, denoted 
by^iS^s), for j — 0,1,2. The procedure is extended in 



Assume S m = {si, ... s^} is a set of sampling points 
on an irregular grid and X*(S nl ) = {A" 1 *, . . . , X^} is the 
respective vector of measurement. On an irregular grid, 
the translation symmetry of the lattice is lost. The con- 
tinuum FGC functional is then a more suitable model. 
For practical purposed, a tractable approximation of the 
continuum model is needed. 

In [lij . approximations for the sample averages of 
5i(s) and S% (s) are formulated in terms of kernel av- 
erages of the data values. In the following we use the 



notation: 



U{d + 2), 



2d(d-l), {Ai d ). 



„(i) 



d, 4 2) = Ad 2 and 



i , . , — — : — : : — , where the 

summation is over both indices = 1, . . . , N) denotes 
the average of the quantity Aij , weighted by the kernel 
Kfc(r) with bandwidth parameter h. The pair distance is 
denoted by s 



|s.j — Sj||, where ||r|| is the Euclidean 
norm of the distance vector r, and the field increment by 
Xtj =X*( Si )-X*(s. 
constraint is given by 



1,3 



X*j = X*(si) — X*(sj). Then, the generalized gradient 



„(i) 



(2) 



where a\ — (s 2 j) h ■ Sensible estimates of the spacing 



ai should account for the grid topology. E.g., let 23o be 
the set of near-neighbor vectors of all the points in S m . 
If Q5o contains No vectors and Aj denote the lengths of 
the vectors in 25q, then af = Yli^i ^t- Similarly, if 



the lattice constant, and the bandwidths hi and h 2 that 
determine the range of influence of the averaging kernel. 
The latter are determined from the consistency principle 
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s 2p \ 



where p — 1, 2. 



= (si tj ) h3 , h 3 = V2h 2 , h 4 = 2h 2 , the generalized tem £x(sj - sj) Xj + fx = C x (si - z ; ), Vi = 1 



SPATIAL PREDICTION 

Let Z p = {zi, . . . zk} be a set of prediction points, dis- 
joint from S m , Vi = S m U {z;}, and V = Z p U S m . The 
predictions will be denoted by {X\(zi), I = 1, ...,K} 
and the respective prediction vector by X\(Z p ). The in- 
crements corresponding to Vi will be denoted by a p (Vi), 
p = 1,2. Typically, single-point prediction is applied 
sequentially over all points in Z p . Multiple-point predic- 
tion is possible in the SSRF framework, but this letter 
focuses on single-point prediction. 



Optimal Linear Prediction 

In geostatistics, spatial prediction is based on the Best 
Linear Unbiased Estimator (BLUE), commonly known 
as Kriging Different variants of kriging exist, de- 

pending on the hypotheses about the normality of the 
data and the behavior of the mean. These methods are 
generalizations of the linear minimum mean square er- 
ror (LMMSE) estimators, also known as Wiener filters 
|l6| . Ordinary kriging (OK) is the most common variety. 
It is applied to normally distributed data, with an un- 
known mean that can be considered as locally constant. 
A single-point prediction is obtained as a superposition 
X(zi) = J2j=i Aj X(sj), where Sj are all the points inside 
a local search neighborhood, B(si), around z/. The pre- 
diction error is defined as e{si) = X(zi) — X{z{). The op- 
timal linear coefficients should minimize the mean square 
error conditional on the zero-bias constraint 53j=i Aj = 

1, i.e., the expression E [e 2 (s;)] +^1 \ J2jf-i\j — , where 
fi is a Lagrange coefficient. This leads to the linear sys- 

,Af 



curvature constraint S2 (s) is given by 



5 2 (s) = 



,(2. 



Mi 



( X tj 



„(3) 



/'2 



,.(1) 



(3) 



where [i\ and (12 are 1 + o(e) constants that depend on 
the sampling network topology. The hi, fi2 are defined so 
as to satisfy asymptotic bias and consistency properties 
[lfij . They introduce explicitly in the problem of model 
inference four parameters linked to the topology of the 
sampling network: the spacings a\ and a 2 that replace 



and YljLi Aj = 1> where Cx(r) is the centered covari- 
ance function. OK is an exact interpolator, meaning 
that {X(si) = A*(si),Vsi e S m }. Exactitude is not al- 
ways desirable, since it ignores measurement errors and 
over-constrains the predictions. For a fixed-size search 
neighborhood, the numerical complexity of an M-point 
OK prediction is 0(K M 3 ). 



Spatial Prediction based on FGC Functional 

Once the parameters of the FGC model have been de- 
termined from the data, prediction of the SRF at z; is 
possible by means of at least two different approaches. 
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First, the corresponding covariance function is deter- 
mined and spatial predictions are then obtained using the 
OK predictor. In this case, the only difference introduced 
by the SSRF functional is the covariance estimator. 

Here we propose a different predictor, obtained 
by maximizing the conditional probability density 
fx[X x (xi)\X x (S m )]. Since fx[X x {xi)\X x (S m )] = 
fx [X x (Vi)}/f x [X\(S m )], the prediction is obtained by 
maximizing fx [X\(Vi)]. In principle, this requires solv- 
ing the equation SH{ gc /5X\(zi) = 0, where gxl(zi) * s * nc 
variational derivative of the functional given by Eq. l|T]) 
with respect to X\(zi). In practice, Hf gc is replaced by a 
discrctized estimator, Hf gc - Since Hf gc is a bilinear func- 
tional, the prediction follows from the solution of the 
linear equation: 



dH [sc [X x (Vi)] 



dX x ( Zl ) 



(4) 



X A ( Zi ) 



Hf gc is obtained by means of the estimators {Sj(s),j — 
0,1,2}, using the ergodic hypothesis JdsSj(s) « 
\ Sj(s), which leads to: 



\n\ 



(5) 

where the spatial averages involve the sampling points 
and the prediction point as well. In light of J2J), and 
(JHJ), equation Q leads to the following linear predictor 



(6) 



where q% = q 2 = 1, = = — 1, (X*) l h is the kernel 
average of the sample, centered at the prediction point 

\ A ih — r ¥ /„ — ' y') 



and the linear coefficients {cj, i = 1, 2, 3, 4} are given by 

Ci(yi) = bi(yi) (N+i) ghl (yi), = c^ m il b 2 (Vi) = 

4> = C/a P (Vi) and 

Ei K ft( s i - Si) 



9h(Vi) 



The summation in gh(Vi) extends over all the N(N— 1) /2 
non-identical and non-repeating pairs of sampling points. 
Defining the linear weights 

iy l> 1 + Cl (Vi) + c 2 {Vi) - c 3 (Vi) - Ci (Vi) ' 
the prediction is expressed as 



xx(z l ) = J2 4 p=1 x p(Vi)(x*y hp 



(8) 



Properties of the FGC Mode Predictor 

The present formulation of the FGC mode predictor 
(FGC-MP) is closer to simple kriging than OK, since the 
mean is assumed to be known. However, unlike simple 
kriging the mean does not have to be constant, provided 
that it changes slowly so that the energy contributions 
due to the square gradient and curvature of the mean can 
be ignored compared to the fluctuations. In this respect, 
the predictor resembles OK, which allows for slow (but 
unknown) variation of the mean. Predictions of the FGC- 
MP are independent of r) , while the prediction variance 
is linearly proportional to 770- 

The FGC-MP is linear and unbiased. Since the joint 
FGC pdf is Gaussian, the mode estimate is equivalent to 
the minimum mean square estimate. Thus, the FGC-MP 
is an optimal linear predictor. The main differences with 
OK result from the use of the energy functional in the 
FGC SSRF: (1) The FGC-MP is not an exact interpo- 
lator, because it does not use the data at the prediction 
point. (2) The single-point FGC-MP provides an explicit 
expression for the prediction, while kriging requires solv- 
ing a linear system. (3) The FGC-MP does not require 
specifying a search neighborhood around the prediction 
point; in kriging definition of a search neighborhood re- 
quires an iterative procedure based on cross-validation of 
the predictions with the data. (4) The FGC-MP incor- 
porates two sets of parameters: the first set determines 
the spatial dependence of the SRF, while the second set 
depends on the topology of the sampling network. The 
influence of the sampling topology is not explicitly ac- 
counted for in kriging. (5) The uncertainty estimate 
involves the SSRF covariance function, for which there 
are no explicit solutions in d — 2, unlike d = 1,3 [lij . 
Obtaining the covariance in d = 2 requires performing 
numerically a univariate (for isotropic dependence) inte- 
gration of the spectral density. (6) Regarding multiple- 
point estimates the FGC-MP has a numerical complexity 
0(K 3 ), derived from solving a linear system of K coupled 
equations at the prediction points, while the numerical 
complexity of kriging is 0(K M 3 ). 



PREDICTION USING SIMULATED SAMPLES 

At 400 randomly distributed points on a square do- 
main of length L = 100 we generate 100 independent 
"samples". These represent realizations of a Gaussian 
random field with m x = 50, and an exponential covari- 
ance function C x (r) = a\ exp(— ||r||/6 e ), where er x = 10, 
and h e = 10. The Cholesky LU decomposition method 
is used for the simulations. We partition the 400 points 
into a training set S m of 100 randomly selected points, 
and a prediction set, Z p , including the remaining points. 
We use the first set to determine the optimal SSRF pa- 
rameters, and then predict the values of the field at the 
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TABLE I: Statistics of OK performance. 





]VI i n i mil m 


]VI ax i mil m 


iviean 


JVIedian 


bias 


-1.90 


2.10 


-0.04 


-0.14 


mae 


5.00 


6.95 


6.03 


6.11 


rmse 


6.28 


8.98 


7.76 


7.81 


mare 


0.10 


0.15 


0.12 


0.12 


rmsre 


0.13 


0.25 


0.17 


0.17 


R 2 


0.38 


0.76 


0.62 


0.63 


TABLE II: Statistics of FGC-Mode performance. 




Minimum 


Maximum 


Mean 


Median 


bias 


-1.58 


2.40 


-0.03 


-0.00 


mae 


5.07 


7.07 


6.13 


6.10 


rmse 


6.40 


9.04 


7.73 


7.73 


mare 


0.10 


0.16 


0.13 


0.13 


rmsre 


0.13 


0.27 


0.18 


0.17 


R 2 


0.37 


0.76 


0.60 


0.61 



locations of the prediction set. The triangular kernel is 
used in the FGC-predictor mode. Predictions are also 
generated using the Ordinary Kriging method. 

The performance of the predictors is evaluated using 
the bias, the mean absolute error (mae), the root mean 
square error (rmse), the mean absolute relative error 
(mare), the root mean square relative error (rmsre) and 
the linear correlation coefficient (i? 2 ). The means are cal- 
culated with respect to the values at the 300 prediction 
points. Statistics of these quantities over the 100 sam- 
ples are shown in Tables [I] and [H] The kriging predictor 
is applied with the a priori parameters of the exponen- 
tial covariance (instead of the inferred covariance model 
from the data). This choice aims at testing the FGC- 
Mode Predictor against the "true" model. The results 
show that the two predictors perform very similarly. 



CONCLUSIONS 

A fast linear optimal predictor, with applications in 
the analysis of spatial data, is proposed. The predictor 
is based on generalized random fields which represent the 
spatial dependence in terms of interactions. An explicit 
expression for single-point prediction is obtained. The re- 
duced numerical complexity of the FGC-Mode predictor 
may promote the use of cross-validation procedures for 
model parameter inference, instead of the commonly used 
parametric methods. The SRF representation, which 
is based by construction on an objective function, pro- 
vides a unified framework for model parameter estima- 
tion, spatial prediction and constrained (respecting the 
data) simulation. This is in contrast with classical ap- 
proaches that require the introduction of ad hoc objective 
functions 0,0 for simulations (e.g., by means of simu- 



lated annealing.) Expressions for the prediction uncer- 
tainty and a linear system for multiple-point prediction 
have also been derived and will be reported elsewhere 
[l7j . The multiple-point predictor accounts for interac- 
tions between the prediction points that may lower the 
total "energy". Such interactions are missed in single- 
point prediction. Finally, the FGC focuses on short-range 
interactions, but long-range dependence can be incorpo- 
rated in the SSRF framework with suitable modifications 
of the energy functional. 



Acknowledgments 

This research is supported by the Marie Curie TOK 
Action of the European Community (project "SPAT- 
STAT" MTKD-CT-2004- 014135) and co-funded by 
the European Social Fund and National Resources 
(EPEAEK II) PYTHAGORAS. 



t 

[1 
[2 
[3. 
[4 
[5 
[6 
[7 
[8 

[9 
[10 

[11 

[12 

[13 
[14 

[is; 

[16 
[1? 



Elect ronic address: |dionisi@mred.tuc.gr | 

URL: http : //www.mred. tuc . gr/home/hristopoulos/dionisi .html 

Electronic address: elogne@mred.tuc.gr 
P. K. Kitanidis, Introduction to Geo statistics: Applica- 
tions to Hydrogeology (Cambridge, 1997). 
Y. Rubin, Applied Stochastic Hydrogeology (Oxford Uni- 
versity Press, New York, 2003). 

M. E. Hohn, Geostastistics and Petroleum Geology 
(Kluwer, Dordrecht, 1999). 

H. Hamzhepour and M. Sahimi, Phys. Rev. E, 73, 056121 
(2006). 

G. Christakos, Random Field Models in Earth Sciences 
(Academic Press, San Diego, 1992). 

P. Goovaerts, Geostatistics for Natural Resources Evalu- 
ation (Oxford, NY, 1997). 

G. Christakos, and D. T. Hristopulos, Spatiotemporal En- 
vironmental Health Modelling (Kluwer, Boston, 1998). 

G. Winkler, Image Analysis, Random Fields and Dy- 
namic Monte Carlo Methods: A Mathematical Introduc- 
tion (Springer, New York, 1995). 

A. Leow et ai, Neurolmage, 24, 910 (2005). 

M. Yaglom, Correlation Theory of Stationary and Related 

Random Functions I (Springer, New York, 1987). 

C. Lantuejoul, Geo statistical Simulation: Models and Al- 
gorithms (Springer, Berlin, 2002). 

H. Wackernagel, Multivariate Geostatistics (Springer, 
Berlin, 2003). 

D. T. Hristopulos, SIAM J. Sci. Comput. 24 2125 (2003). 
D. T. Hristopulos and S. N. E logne, submitted to IEEE 
Trans. Inform. Theor., cs.IT/0605073 

S. N. Elogne and D. T. Hristopulos (2006). 
math. ST/0603430 

J. Ruiz-Alzola, C. Alberola-Lopez and C.-F. Westin, Sig- 
nal Processing, 85(2), 413 (2005). 
D. T. Hristopulos and S. N. Elogne, in preparation. 



